Monte Carlo Sim for Power.R

In jimmyrigby94/rrs: rrs: Robust Response Surface Analysis

library(rrs)
library(tidyverse)

# Defining Power Assumptions ----------------------------------------------

# Define cor_mat
cor_mat<-matrix(c(1, 0,
                  0, 1), 
                byrow = TRUE, 2,2)


# Defining betas
beta<-list(c(0, 0, -.075, -.075, .15),
           c(0, 0, -.05, -.05, .1),
           c(0, 0, -.025, -.025, .05))

# How many iterations should be run per sample size?
iter<-1000

# Increments of sample size ranges
sample_range <- seq(100, 1000, by = 50)

# Initalizes Results
power_results<-data.frame()

# iterates through betas
for(b in 1:length(beta)){
  
  message("Beta: ", b, " of ", length(beta))
  
  # Finds appropriate error for each beta
  error_sig<-find_sig(n = 1000, cor_mat = cor_mat, beta = beta[[b]], target_var_y = 1)
  
  # iterates through sample sizes
  for(s in 1:length(sample_range)){
    
    message("Sample Size: ", s, " of ", length(sample_range))
    # Initializes iteration data frame starting new each time
    output<-data.frame()
    
    # Runs iter iterations Monte Carlos samples for each beta sample size combo
    for(i in 1:iter){
      
      # simulates data for iteration i
      sim_data<- gen_response_surf_x(sample_range[s], cor_mat = cor_mat, x_names = c("precep_help", "nurse_help"))%>%
        gen_response_surf_y(beta = beta[[b]], sigma = error_sig, y_name = "mastery")
      
      # Runs response surface model
      m<-resp_surf(dep_var = "mastery", fit_var = c("precep_help", "nurse_help"), data = sim_data, robust = FALSE)
      
      # Pulls out lines of interest
      output<-bind_rows(output, m$loi)
    }
    

     # Calculates Power Results for model
     power_results<- bind_rows(power_results, output%>%
                                 group_by(Line_of_Interest, Parameter)%>%
                                 summarise(Estimate = mean(Estimate), 
                                           Power = mean(P_value<.05))%>%
                                 mutate(`Sample Size` = sample_range[s],
                                        `Effects` = paste(beta[[b]], collapse = "_"))
     )
     
  }
}

# Verifying the approriateness of the effect size by calculating R-squared.
sigs<-c()
  for(b in 1:length(beta)){
    sigs[b]<-find_sig(n = 1000, cor_mat = cor_mat, beta = beta[[b]], target_var_y = 1)
  }

1-sigs


p<-power_results%>%
  filter(Parameter == "Quadratic", Line_of_Interest == "Line of Incongruence")%>%
  mutate(Effect = case_when(Effects == "0_0_-0.025_-0.025_0.05" ~ "Small",
                            Effects == "0_0_-0.05_-0.05_0.1" ~ "Medium",
                            Effects == "0_0_-0.075_-0.075_0.15" ~ "Large"))%>%
  ggplot(aes(x = `Sample Size`, y  = Power, lty = Effect))+
  geom_smooth(se = FALSE, color = "black")+
  theme(panel.grid = element_blank(),
        panel.background = element_blank(),
        axis.line.x = element_line(),
        axis.line.y = element_line())+
  geom_hline(yintercept = .8, lty = 4)+
  geom_vline(xintercept = 202, lty = 1)+
  geom_vline(xintercept = 407, lty = 2)+
  labs(title = "Estimated Power Curves for Quadratic Effect Along Line of Incongruence",
       subtitle = "Generated using 10,000 Monte Carlo Samples per Condition",
       caption = "Note. Sample size was changed by increments of 50. Curves are generated by a loess smoother.")+
  scale_x_continuous(breaks = seq(100, 1000, by = 50))

p

sim_help<-gen_response_surf_x(1000, cor_mat, x_names = c("precep_help", "nurse_help"))%>%
            gen_response_surf_y(beta = beta[[1]], y_name = "mastery")

m<-resp_surf(dep_var = "mastery", fit_var = c("precep_help", "nurse_help"), data = sim_help, robust = FALSE)

plot_ly_surf(m, 2, -2, 2, -2, inc = .25, xlab = "Protege Help", ylab = "Mentor Help", zlab = "Task Mastery", showscale = TRUE)